Numerical modeling of collisional dynamics of Sr in an optical dipole trap 
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We describe a model of inelastic and elastic collisional dynamics of atoms in an optical dipole 
trap that utilizes numerical evaluation of statistical mechanical quantities and numerical solution 
of equations for the evolution of number and temperature of trapped atoms. It can be used for 
traps that possess little spatial symmetry and when the ratio of trap depth to sample temperature is 
relatively small. We check parameters calculated numerically against analytic calculations for power- 
law traps and compare results with experiments on **Sr, which has well-characterized collisional 
properties. 



I. INTRODUCTION 

Understanding the collisional dynamics of trapped, ul- 
tracold atoms is essential for optimizing forced evapo- 
rative cooling [l|, and obtaining quantum degenerate 
Bose and Fermi gases It also allows determina- 
tion of ultracold collision properties from the evolution of 
number and teniperature in a trapped sample of atoms 
or molecules 0, [a, 0| • 

Many recipes have been presented for relating the evo- 
lution of the trapped gas to underlying physical parame- 
ters. Typically the collisional dynamics are described by 
differential equations for the time rate of change of the 
atom number N and total energy E, as originally sug- 
gested by [Hi]. The method has been extended and de- 
veloped in many other works [1, S U^, [H [H, [H, Q [ll ■ 
The standard treatment of evaporation is described by 
Luiten et al. [To| . which derives expressions for ther- 
modynamic quantities from the kinetic equations using 
an assumption of sufficient ergodicity and a truncated 
Boltzmann velocity distribution. Analytic evaluation of 
these expressions is straightforward for power-law traps. 
Noteworthy subsequent improvements over this work in- 
clude the addition of effects of time-dependent potentials 
PH . ener gy-d ependent cross sections [l^, and quantum 
statistics Prescriptions have been offered for opti- 

mizing evaporation |14l | and deriving scaling laws [l5|. 
Direct Monte Carlo simulations have also been presented 
to relax the assumption of sufficient ergodicity [T6| and 
treat hydrodynamic effects [17]. 

A common simplifying assumption is that 77, the ra- 
tio of trap depth et to sample temperature fc^T with 
Boltzmann constant fcs, is large. For example, this 
yields analytic expressions for thermodynamic quantities, 
and allows approximation of optical dipole traps Tsl] as 
parabolic potentials [l^ . By taking advantage of the high 
degree of spatial symmetry in a linear potential, analytic 
expressions for thermodynamic quantities were derived 
for the low-?7 situation (77 < 4) in this particular geome- 
try [2^. It is worth emphasizing that Luiten's model [13] 
is, in principle, valid for low 77, as long as the assumptions 
of ergodicity and a truncated Boltzmann distribution are 
also valid. 

If the potential lacks the ideal shape of a power-law 



trap, simple analytic expressions for many quantities of 
interest cannot be found, and numerical methods are re- 
quired. This is the case for low-77 in an optical dipole 
trap and especially when gravity is significant. *^Sr in 
an optical dipole trap falls into this situation because 
of its large mass and extremely small s-wave scatter- 
ing length a — — 1.4(6) ao [IH, where the Bohr radius 
flo = 5.3 X 10~^^m. In the present work, we present nu- 
merical methods for describing collisional dynamics in a 
static trap in this regime, based on the works of Luiten 
et al. |lCl| ] and Comparat et al. [l^, and apply them to 
describe experiments with ®*Sr. 

This paper is organized as follows. Section ITIl describes 
the experimental setup, and then section Hill presents the 
collisional processes important in the trapped sample and 
the differential equations for evolution of N and E. As a 
check of numerics, and to study ranges of validity of var- 
ious approximations, a comparison of numerical calcula- 
tions with corresponding analytic results for power-law 
traps is presented in section HVl A brief introduction to 
the numerical code, as applied to trapped ^^Sr, is pro- 
vided in section |Vl Finally, the elastic collision cross sec- 
tion of '^^Sr inferred from the fit of the numerical model 
to observations is compared with theoretical results in 
section I VII 



II. EXPERIMENTAL SETUP 

The creation of samples of ^*Sr atoms in an optical 
dipole trap (ODT) starts with laser cooling and trap- 
pin g phases that have been described in detail previously 
p^l23l[2l]. Atoms are trapped in a magneto-optical trap 
(MOT) operating on the 461 nm ^Sq-^Pi transition (Fig. 
[Ij and cooled to about 2mK. There is a decay channel 
from the ^Pi state, to the state with a branching 
ratio of 2 X 10~^. This can subsequently populate the 
metastable level. atoms are repumped by ap- 
plying a 3/im laser resonant with the ^P2-^D2 transition 
that returns these atoms to the ground state. The re- 
pumped sample of atoms contains about 3.5 x 10* atoms. 

After this initial MOT stage, the 461 nm light is ex- 
tinguished and the atom sample is transferred with more 
than 50% efficiency to a second MOT operating on the 
^So-^Pi intercombination line [1^. The sample is cooled 
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FIG. 1: (Color online) Atomic Sr energy levels involved in the 
two-photon PAS experiments. Decay rates (s~^) and excita- 
tion wavelengths are given for selected transitions. Laser light 
used for the experiment is indicated by solid lines. Atoms de- 
caying to the level may be repumped by 3 jj.m light. 



to producing densities of 10^^ cm 

Atoms are then transferred to an ODT generated from 
a 21 W, 1064 nm, linearly-polarized, single-transverse- 
mode fiber laser. The experimental setup in this work is 
shown in Fig. [21 The trap is in a crossed-beam configura- 
tion, derived from the first order deflection of an acousto- 
optic modulator. The beam is focused on the atoms with 
a minimum e^^ intensity-radius ofw^ 100 /im. It is then 
reflected back through the chamber to intersect the first 
beam at 90 degrees and refocused to have approximately 
the same waist at the atoms. Both beams lie in a plane 
that is inclined 10.5° from horizontal. 

The number of atoms and sample temperature are de- 
termined with time-of-flight absorption imaging using the 
^Sq-^Pi transition. The ODT trapping potential is cal- 
culated from measured laser beam parameters and the 
polarizability of the ^5*0 state [1^, and it is checked by 
measuring the trap oscillation frequencies through the 
parametric resonance technique [27|. This allows us to 
infer the sample density profile from the temperature and 
number of trapped atoms. 



III. MODEL OF COLLISIONAL DYNAMICS 



The evolution of N and E is described by a system of 
differential equations. Different terms in the equations 
represent physical processes such as elastic and inelastic 
collisions, and processes involving laser fields. 




FIG. 2: (Color Online) Schematic of our experiment illustrat- 
ing the overlap of ODT beams with MOT beams and relative 
positions of magnetic coils. 

A. Description of Basic Processes 

1. Background collisions and inelastic collisional losses 

One-body losses due to collisions with background gas 
and two- and three-body inelastic collisional losses are 
described by the local equation 

h = -TbgU - /3mn^ - Ln^, (1) 

where n is the atomic density. In simulations described 
here, we will assume that the loss rate constants Tbg, Pim 
and L are independent of temperature. Integrating Eq. 
[T] over the trap volume gives 




where the effective volumes are 

Vg ^ J— [ d^nirW. (3) 

peak 

We have also made use of the relationship between peak 
density and total number, npg^akVi = 

The energy or temperature evolution due to these pro- 
cesses for a constant trap potential can be found as fol- 
lows. The rate of energy change in an infinitesimal vol- 
ume dV is 

dE = ~h{r)dV{U{r) + ^^(r)), (4) 

where U (r) is the trap potential and Ek (r) is the average 
kinetic energy per atom located at r. U (r) is defined to 
have a value of ?7 = at the trap minimum. We assume 
a truncated Boltzmann phase-space distribution (this is 
valid if the trap is sufficiently ergodic [iq|), which implies 
that the kinetic energy in a given differential volume also 
obeys a truncated Boltzmann distribution truncated at 
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the kinetic energy required for an atom to escape the 
trap from the differential volume. Thus the position- 
dependent average kinetic energy can be expressed as 

Integrating Eq.d] yields the rate of change of total en- 
ergy 



We have introduced the effective kinetic energies 
1 



T, 



q — n 

n , 

peak 



d-'r[n{v)]'^Ek{r), 



and effective potential energies 
1 



P, 



q — q 

''^peak 



d^r[n(r)]«t/(r). 



(6) 



(7) 



(8) 



Note that the total energy of atoms inside the trap is 
E = {Ti+ Pi)^. In the high-?? limit, ^fc(r) = ffcsT, 
which can be taken out of the integrals to yield 

V2 9 V3 



E 



T,,^^N + p.,,,^N^ + L^N^y (9) 



2. Off-resonant laser scattering 

The scattering of off-resonant photons, such as from 
the ODT laser, heats the atoms due to momentum dif- 
fusion (MD) [2^ . The rate of change of total energy due 
to this process is 



E 



MD 



— ^laserNE. 



recoil 5 



(10) 



where the total scattering rate of light from the laser 
field is given by Tiaaer- If the light scattering rate is 
dominated by one transition, Via, 



_ S07/2 



which So is the saturation parameter, 5 is the detuning 
of the laser, and 7 is the linewidth of the transition. The 
recoil energy is ErecoU = kBTrecoU = h^k^/m where h 
is Planck's constant h divided by 2tt, k is the photon 
circular wavenumber, and m is the atom mass. 



3. Evaporation 

To describe the rate of atom loss due to evaporation, 
we follow the treatment of [lo| , which assumes ergodicity 



and a truncated Boltzmann distribution, 

'^0 „..„r /rr/„N , 2 



/(r,p) 



(27rmfcBT)3/2 



exph([/(r)+pV2m)/fcBr] 

(11) 



where no is the peak density (the density at the trap min- 
imum) in an infinitely deep trap, &{x) is the Heaviside 
step function, and p is the atom momentum. This yields 
a density distribution given by 



n(r) = npeafeAe-^W/^«^ierf 




^ et-Ujr) 

where the normalization constant A is given by 
no 



, (12) 



|erf 









-2 



exp 



keT 



(13) 

It is important to note that no is no longer the peak 
density since the peak density is redefined as 



ripeak = "-(r)|c/(r)=0 

~ no< erf 



^JtilksT 



-2\J et/TrfesTexp 



- et/ksT 



(14) 



The total number of atoms lost per unit time due to 
evaporation can then be written as 



Nev = -TevN, 



where the evaporation rate per atom is 
N 



(15) 



(16) 



Here, dei is the elastic collision cross section, which is 
assumed to be collision-energy independent in this treat- 
ment. V — (^f^)^^^ is the mean atomic velocity, and 
the effective volume for elastic collisions leading to evap- 
oration is 



keT 



where 



/ * depie)[{et - e - kBT)e-'/''^^ + keTe-^ 
Jo 

(17) 



A = (27r?iV™fcsT)i/2 



(18) 



is the thermal de Broglie wavelength. The density of 
states in the trap is given by 



P(e) = 



27r(2rn)3/2 



U{v)<et 



d^r^Je~U{v). (19) 
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Similarly, the rate of change of total energy due to 
evaporation is 

Eev — —^evNEev, (20) 

where the average energy loss per evaporated atom is 



Eev — 



(21) 



with 



A3 

kf,T 



I ' dep(e) [fc^Te-^/^^^^ - (e* - e + kBT)e-'^] . 
Jo 

(22) 



4- Final equations 

Accounting for all processes, the equations for number 
and energy evolution become 



N = 



E 



(23) 



Vi 



-L ( "^^3 ) + T laser NErecoil 

'—^Aaeive 'Vev e* H — kbi 



Ve, 



.(24) 



B. Ratio of Elastic to Inelastic Cross Sections 

For interpretation of some experimental data, it is nec- 
essary to calculate additional factors that describe the 
roles of inelastic and elastic two-body collisions. As 
pointed out in [20], losses due to inelastic two-body col- 
lisions and evaporation both scale with the number or 
density of atoms squared, which makes it difhcult to mea- 
sure the rates of the two processes individually from loss 
curves alone. Additional measurements of temperature, 
along with a statistical mechanical model of evaporation 
and inelastic loss in the trap can overcome this problem. 

When all loss processes except evaporation and two- 
body inelastic collisions are negligible, because inelastic 
loss is a heating mechanism and evaporation is a cool- 
ing mechanism, the atoms tend toward a steady state 
temperature that is independent of the number of atoms 
in the trap [13] ■ The ratio of trap depth to the steady 
state temperature is the equilibrium rj (at constant trap 
depth). Assuming that the steady state temperature is 
reached, the local loss equation is n = — /3n^ where (3 is 
the total two-body loss rate constant (including inelastic 
and evaporative losses) (5 = Pin + f{r])(3ei, where /(?]) 
is the fraction of an atom's elastic collisions that result 



in the evaporation of that atom [3l|, and (3in/ei ^-re the 
inelastic and elastic collision rate constants. For an un- 
truncated Boltzmann (or very large rj) j3ei — \plaeiv, 

where v 



Typically, A„ = 2K, where K is the 
inelastic collision event rate constant, each collision re- 
sults in the loss of both atoms, and — \/2(TinV. The 
relative contributions of inelastic and elastic losses can 
be described by [20j 



Pin — 
Pel = 



1 



/7+1 

7 



/7+1 

where the ratio of elastic to inelastic cross sections is 



(25) 
(26) 



7 ^ 



Pel 
0in 



<7el 



E-Ed 
f{Eev - E) 



(27) 



and E = (ri+Pi)/^! is the average energy of the trapped 
atoms, Ed = {T2 + P2)/V2 is the average energy of atoms 
lost due to inelastic collisions, and Eev is the average 
energy of atoms lost due to evaporation (Eq. [H]). All 
variables on the right-hand side depend only on rj, which 
is constant at equilibrium. 

The calculation of / is somewhat involved. It is given 
in terms of other coUisional parameters by 



(28) 



where Fei, is the evaporation rate per atom (Eq. I16p and 
Tel is the elastic collision rate per atom (number of col- 
lisions per second, on average, that an atom in the trap 
suffers). Several ancillary quantities must be calculated 
in order to find Tei ■ The number of elastic collisions per 
second per unit volume at position r is [l^, [s^] 

^(r) = ^ J Spi j d^P2cre;lvi - V2l/(r,pi,i)/(r,P2,t), 

(29) 

where Vi = pi/m, V2 — P2/?7T-, and we have assumed an 
energy-independent elastic cross section. For a truncated 
Boltzmann distribution independent of time, this integral 
must be evaluated numerically. In the large-?; limit, it 
reduces to the standard result 



Z{U{v))^2aei[n{v)\' 



ksT 



(30) 



The total number of elastic collisions happening per 
second inside the trap can be found by integrating Z over 
the volume of the trap. There are two atoms involved in 
every collision. So the average rate at which an atom 
undergoes elastic collisions in the trap is 



Tel = 



N 



(31) 



and fTeiN is the number of evaporating atoms per sec- 
ond. For an untruncated Boltzmann distribution, Tei in a 
square, a parabolic, and a linear potential are ^/2no(Jeiv, 
noaeiv/2, and noaeiv / 4:\/2 respectively [20j ]. 
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IV. CHECK OF NUMERICAL QUANTITIES 



Effective Volumes V„ 



In principle, all the quantities necessary to solve Eqs. 
[23l and [24] and model coUisional dynamics can be calcu- 
lated, but in practice approximations are necessary to ar- 
rive at analytic results. This is straightforward for high-ry 
conditions [y], and also in situations of low- 77 and with 
sufficient trap symmetry pol |. For low- 77 conditions and 
traps that lack spatial symmetry, numerical evaluation of 
statistical mechanical quantities is the only option, and 
that is the approach we follow. 

To provide a check of our numerical evaluations, we 
compare numerical results (numeric, general (NG)) with 
analytic forms for power-law traps. For the NG spatial 
integrals for p{e), Vq, Tq, and Pq, the integration extends 
over the region in the trap determined from an algo- 
rithm that finds the saddle points of the potential. An 
interpolating function representing p{e) is used in evalu- 
ation of Vev and Xpv in integrals over an energy interval 
from to et . An adaptive-step-size integration routine in 
Mathematica is used to evaluate these integrals. 

We use two different analytic forms, one in which the 
standard high-77 approximations are used (analytic, high- 
77 (AH-?])) [9|,[20| and one in which the high-?] approxima- 
tion is not made and a truncated Boltzmann distribution 
is assumed (analytic, general (AG)) p^ . 



A. Power-Law Traps 

Ideal power-law (PL) trap potentials have the form 



U{r) = et 



Rr. 



(32) 



For the linear and parabolic traps j — 1 and j — 2 
respectively. The square potential can be viewed as a 
"f —f 00 limit 



C/(r) 



Square 



cx), for r > Rq 
0, for < r < Ro 



Here, and Rq can be interpreted as scaling parameters. 
Truncated traps, which lead to truncated Boltzmann dis- 
tributions, can be expressed as 



Effective volumes Vq are defined in Eq. [3l From the 
expressions for n{r) and npeak given in Eqs. fT2lfT4l the 
AG form of Eq. [3]can be evaluated exactly as: 



47r 



[P(?7,3/2)]^ 



-Ro 



^2^-qU(r)/kBT 



P 77- 



U{r) 3 
fc^' 2 



with the incomplete Gamma function 



P{-q,a) 



(36) 



(37) 



The denominator in Eq. [37]is the Gamma function T{a). 

and 



Transforming to the unitless variables f 
ill 



EM — •qf'^ yields 



Jo 

[P(?7-?7fT,3/2)]9 



[P(?7,3/2)]^ 



dr. 



(38) 



If we assume large 77, the incomplete gamma function P 
tends to unity, and we arrive at the AH-r; approximation 
for Eq. EH] 



AnRi 



47ri?gr(3/7+l)-(3/7)r(3/7,g77) 

3 (ot)^/'^ 
47rp3r(3/7+l) 



3/7 



r(3/7 + i), 



(39) 



where the approximation step takes the limit of the upper 
(33) incomplete gamma function r(s, 77) = f'^^e^^dt ~> 
when 77 is large. 

Fig. [3] shows comparison plots of Vi, scaled by Vq — 
47rPQ/3. The NG curves are identical to the AG ones as 
expected. The high-77 niodel is a reasonable approxima- 
tion for 77 > 4 in agreement with ^] . 



and 



(34) 



C. Effective Volume for Elastic Collisions Leading 
to Evaporation Vev 

Effective volumes for elastic collisions leading to evapo- 
ration, Vev, for PL traps in the AG treatment are defined 
in Eq. 42a in p| 



where U{Ro) = et is the trap depth and Rq is the bound- 
ary. 



yAG _ A 3 

* ev 



A^Coo[r7P(r?,3/2-|-3/7) 
-(5/2 + 3/7)P(r7,5/2 + 3/7)], (40) 
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FIG. 3: Plots of Vi vs rj for (A) Linear, (B) Parabolic, and (C) Square potentials. 




FIG. 4: Plots of Vev vs rj for (A) Linear, (B) Parabolic, and (C) Square potentials. 



where the thermal de Broghe wavelength A is defined in 
Eq. [TH and 

Coc = ApLr(3/2 + 3/7)(fcBr)3/2+3/7. (41) 

for PL traps can be expressed by 



A 



PL 



3/2 



Vo r(3/7 + 1) 



(42) 



Accordingly, Eq. [JD] can be rewritten as 

Vei'^h) = Forr'/^r(3/7+l)[77P(r/,3/2 + 3/7) 

-(5/2 + 3/7)P(?7,5/2 + 3/7)]. (43) 

Because P tends to unity at high-?7, the AH-?7 approxi- 
mation for Eq. [33] is 

Vei"-'''{l) = VoV-'/mVl + 1)N - (5/2 + 3/7)]. (44) 

Comparisons plots of Vev, scaled by VoT]^^^'^ are pre- 
sented in Fig. m The NG and AG curves are identical as 
expected. 

The volume parameter X^v (Eq. for PL traps in 
the AG model is given by 



(45) 



which can be written as 

XfJ'h) = Fory-3/7r(3/7 + 1)P(^, 7/2 + 3/7). (46) 



When r] is high, this becomes 



X. 



AH-r] 



(7) = T/o?7-3/7r(3/7+l). 



(47) 



Comparisons of X^y, scaled by Vb?/^^^^ are plotted in 
Fig. [51 The agreement between the NG and AG curves is 
excellent. However, in Fig. [5] (A), (B) and (C), the AH- 77 
curves match the AG ones for 77 > 10, 77 > 8 and 77 ^ 6 
respectively. 



E. Evaporation Rate Per Atom Fev 



The evaporation rate per atom T^v is defined in Eq. \W\ 



N 



770 



Tev = -j^A'^a^ive "^V^v = '^Au^ive "^Vev (48) 
Eq. [351 in the AG model can be evaluated exactly as 



et) — yAG ^'^elVe , 



(49) 



with V^'-^ and V^'-^ given by Eqs . [38l and [43l respectively. 

In the high-7; approximation, the normalization con- 
stant A and the incomplete gamma function P both tend 
to unity. Substituting V^"^"'"^ and K^^"** from Eq. ^ 
and Eq. [44| into Eq. [48| yields 



ev 

while [§| gives 



YAH-ti 



5 3 

5 + 7 



(50) 



(51) 
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FIG. 5: Plots of Xev vs rj for (A) Linear, (B) Parabolic, and (C) Square potentials. 
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Thus we get the ratio of T^^f^ ^ given in Eq. [50] to large rj. 
those in Eq. (HUfor PL trap potentials: 



(r/ — 11 / 2) /rj, hnear; 

(77 — 4)/ T], paraboHc; 
(77 — 5/2)/?7, square. 



(52) 



Figure [6] compares Tev, scaled by r]e~'^noaeiv in differ- 
ent models. The agreement between NG and AG curves 
is also satisfying. However, it is obvious that, in Fig. [5] 
(A) and (B), the AH- 77 curves match the AG ones for 
77 ^ 8 and 77 > 6 respectively. It is worth emphasizing 
that the ratios of F^^"'' in Eq.[52]only tend to 1 for very 



m 

Tel 



F. Elastic Collision Rate Per Atom Fei 

The elastic collision rate per atom T^i is defined in Eq. 

2A 



2 

N 



(i^rZ(C/(r)) 



d-^rZ([/(r)), (53) 



where Z{U{r)) from Eq. 
rewritten as 



1^1 "0 JVo 

in the AG model can be 



Z(C/(r)) = 



477^ l2haT ,^1 r i 

^J^^CTel dp J dp2e-^P^+^in derSmerJpl+pl-2p,p2COSt 

y m Jo Jo Jo V 



xplpje' 



Pi 



knT 



Pi 



(54) 



with the unitless variables pi = pxj \J2mkBT and p2 = 
P2/ y/2mkBT . 

Substituting Z{U{y)) from Eq . [5^ and simplifying the 
integral with respect to Or, we can rewrite Eg. [551 



rf.^(7) = 



?>2Rlnoa,.iA [2k^ ^^^2^-2^,^ 



dpi 







,~ P1P2 

dp2- 



Jo 3 

:[(75i+p2)3-|pi-p2|3]e-(^?+J^^). 



(55) 



addition, A can be set to unity. Therefore, we have 



el 



(7) 



32R'onoaei 2ksT r ^--2^-2,,'' 



X / dpi / #2^-^ 



x[(j5i+p2)-^-|pi-p2|']e-(P'+^^) 



^^1-'''"'' (7) 7(2^7)=^/^ 
V2 



aeivno 



23/7 



<7eivno, 



(56) 



We note that the integrand in Eq. [55] is suppressed by 
g-2j7r^ and e~^P^~^P^\ so the contribution from small f, 
Pi, and p2 dominates the integral. Accordingly, in the 
high-7; form, the integral can be approximated by ex- 
tending the upper bound of its interval to infinity. In 



which is consistent with the high-7; expressions given in 
Q. Here the last equality is obtained by substituting 
V^"-"^ given by Eq. M 

The NG integration for Fg/ is the most time-intensive 
of all the calculations. To reduce numerical integration 
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FIG. 6: Plots of Tev vs rj for (A) Linear, (B) Parabolic, and (C) Square potentials 
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FIG. 7: Plots of Tel vs 77 for (A) Linear, (B) Parabolic, and (C) Square potentials with different energy resolution in the 
numerical integration. 



time required to make a lookup table for Eq. 1531 we use 



G. Evaporation Fraction f 



PpotenHaliU) = J £r5{U - U {v)) (57) 

to allows us to evaluate Z as a function of U, instead of 
r as 



d\Z{U{T)) = / dUpp,t,r.Ual{U)Z{U). 



(58) 



Then the numerical evaluation of the elastic coUision rate 
is performed as 



Te; = j dUppotential{U)Z{U). 



(59) 



Vei is calculated using a simple constant-step-size rect- 
angle integration method. Interpolating functions for 
Z{U,T) and Ppotentiai(U) are found by numerically eval- 
uating Eqs. [54l and [57l at discrete sets of U and T points 
and U points respectively running from to €4. (The T 
dependence in Z is highlighted here for clarity.) 

The resolution of the grid of evaluation points used to 
form the interpolation of Z was critical. Fig. [7] shows 
the results for different numbers of points evaluated on 
each axis (T and U) N = 2b (dotted hues), 100 (dashed 
lines) and 200 (black solid lines) for comparison. The im- 
provement is considerable with more evaluations. How- 
ever, the computation time varies with iV^, and it takes 
about 24 hours to create the Z interpolating function 
with N = 200, which is deemed sufficient given the bal- 
ance between time and accuracy. 



The evaporation fraction / is defined in Eq. 

fiv) 

For PL trap potentials, reference [9[ gives 



r 

^ ev 



(60) 



{4-\/2?7e linear; 
2rje-'^, parabolic; 
rje^^ / \/2 , square , 



(61) 



which is accurate if t] is very high. Improved expressions 
can be found for high-77 situations following the treatment 
of fH: 



/ 



2\/2[6+2e"(?)-3)+rK4+?))] 
-l+e2'J-2r((l+r)) ' 
AH-r] _ ) e"727r(2»)-3)crf[^]+6v/2») 

e^-J \/27rerf [V2^ -4v^ 

r?e-Vy2, 



linear; 
, parabolic; 
square. 



(62) 



Equation 1621 will tend to Eq. [ni]with very high "q. 

In our formalism, which follows the treatment of [lo| . 
for the general case, Eq. [SOI becomes 



/ 



AG 



-■AG 



rAG ' 
^ el 



(63) 



where F^^*^ and F;^*^ are from Eqs. |49] and [55] respec- 
tively. Following our convention, we obtain a third form 
for the high-?/ approximation by using the high-r/ limits 
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FIG. 8: Plots of f vs rj for (A) Linear, (B) Parabolic, and (C) Square potentials. 



for Vf^~^ and F^j given by Eqs. [50] and [56] respec- 
tively. We arrive at the limiting case that essentially 
follows from [lOj : 

yAH-r) 
— ev 

^ el 




4\/2(r7 - 4)e"'', linear; 
2(7]- ll/2)e-'', parabolic; (64) 
[t] — 5/2)e^''/\/2, square. 



which will also approach Eq. [61] for very high-ry. 

Scaled plots of / given by different methods are shown 
in Fig. [5] The discrepancies of / between AG (circles) 
and NG results (black solid lines) at high-jy correspond to 
discrepancies in T^i due to the finite resolution of the nu- 
merical method as discussed in section llVFI The discrep- 
ancies can be reduced, at the cost of time, by increasing 
the number of sample points when creating the interpo- 
lating function of Z{U,T). Moreover, although all AH-?] 
models give the same results at very high-ry, care must 
be taken using approximation expression for 77 ^ 12. 

V. NUMERICAL SOLUTION FOR A 
REALISTIC TRAP CONFIGURATION 

The great advantage of the numerical code is that it 
can be used for arbitrary trap shapes. This is necessary 
for our experiments with ^^Sr in an ODT because of the 
importance of gravity and the small inclination of our 
trap lasers away from horizontal. For non- ideal traps, it 
is essential to find an expression for the potential, U{r), 
that describes the potential accurately. We employ an 
algorithm to find the trap minimum, the trap depth, and 
the saddle points, which define the trap boundaries and 
allow us to offset the trap so that the minimum is J7 = 0. 
An example is shown in Fig. [9] 

This U (r) is then passed to the numerical integration 
routines for statistical mechanical quantities described in 
Sec. lIIII that are evaluated at a grid of temperature points 
to construct lookup tables. The lookup tables are used to 
find interpolating functions for all quantities, which can 
then be used in place of time-intensive integral evalua- 
tions. Using these lookup tables, given initial conditions. 



the atom number and temperature evolution can easily 
be found from Eqs. and [Ml using an ordinary differen- 
tial equations (ODE) solver in Mathematica. Typically, 
one day is needed to create all lookup tables for a given 
trap potential. After this preparation, a complete ODE 
solution for several seconds of sample evolution only re- 
quires a few seconds of computer evaluation time. 




FIG. 9: (Color online) (A) The x-y cut and (B) the y-z cut 
of the ODT potential for approximately 10 W per beam. The 
boundary of the trap is the set of points with potential en- 
ergy equal to the lowest saddle point, which is along the z axis. 
Gravity is oriented in the -y direction. The beam nearly par- 
allel to the x axis is slightly weaker and less focused than the 
beam along z, leading to the observed asymmetry. 

As described, this procedure is appropriate for a fixed 
trapping potential. It is possible to form lookup tables 
that include the variation of ODT laser intensity, or trap 
depth, and this will allow for changing the potential and 
modeling of forced evaporation. We are currently work- 
ing on this improvement. 

Figure [TU] shows the number and temperature of atoms 
as functions of time for *^Sr ^Sq atoms in the ODT. Var- 
ious quantities such as the one-body loss rate, two-body 
inelastic collision rate constant, and inelastic cross sec- 
tion can be determined by fitting the calculated evolution 
curves to experimental data. 

We exclude the first few hundred milliseconds from the 
fit because we expect that atoms are far from equilib- 
rium and significant population is still trapped in the 
individual beams of the ODT and not in the crossed 
region at this time. For this fit, we obtain the elas- 
tic scattering cross section aei = 260 Oq, one-body loss 
rate Fi = 0.3 s~^, and the ODT photon scattering rate 
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FIG. 10: Variation of atom temperature and number with 
time for **Sr atoms in the ^5*0. The trap used is described in 
Fig.H 



Tout — 0.4 s^^. A good fit is found with = as 
expected since there are essentially no two-body loss pro- 
cesses in this system. 



VI. ENERGY DEPENDENCE OF THE **SR 
ELASTIC COLLISION CROSS SECTION 

Evaporation with an energy-dependent elastic collision 
cross section has been modeled in [l^, but these 
treatments assume a large cross section that varies be- 
cause of the unitarity limit. In **Sr, the cross section 
varies because the scattering length is very small, as 
shown in Fig. [11] [2l[. This variation is significant at 
microkelvin energies, which complicates comparison of 
theory and experiment because a distribution of collision 
energies contributes in a thermal sample in the ODT. A 
full energy-dependent kinetic calculation is beyond the 
scope of this model. We treat the variation in approxi- 
mate fashion by assuming the system can be described by 
an effective, temperature-dependent cross section, {cFei}, 
that is an average of the collision-energy dependent cross 
section. 

To relate (aei) to the underlying energy-dependent 
<yei{Ecoii), first consider the standard expression for total 
collision rate per unit volume in the trapped gas [1^ [s^l 

^(r)==^ j d^pi J d^P2CTe;|vi - V2|/(r,pl,^)/(r,p2,^)• 
(65) 

For simplicity, we use untruncated distributions and as- 
sume a constant density, which yields 



Z = 



27mQ 



dpp^ 



<Tel{Ecoll) 



-e ^MfcBT^ (66) 



Here aei{Ecou) depends on the collision energy EcoU — 
for p = |vi — V2I and the reduced mass fi = m/2. 
For an energy-independent cross section, this reduces to 
[291 



Z = 2(7 elflQ 



(67) 



We assume that the evaporation process is sensitive to 
an average cross section in which the weighting is taken 
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FIG. 11: Measurements and theories for the **Sr elastic colli- 
sion cross section, a^i . The dashed curve is the elastic collision 
cross section 21] for a collision energy of Ecoii = 'iksT, which 
is an average collision energy as described in the text. The 
solid line is an average cross section, weighted by the rate of 
collisions of a given energy, for the temperature T. 



as the contribution of each collision energy to the to- 
tal number of collisions per time in the sample. Using 
this weighting, the average collision energy for a given 
temperature T is 2kBT, and for a given equilibrium tem- 
perature, the average cross section is given by 



Z 



2nl 



(68) 



where Z is calculated using Eg. [66] and the theoretical en- 
ergy dependence of aei [2l| . Contribution to Z is not ex- 
actly equivalent to contribution to the evaporation rate, 
but this is a reasonable approximation in the spirit of 

Fig. [TT] shows variation of (Tei {Ecoii) and (aei) for *^Sr 
as well as experimental data in which {(Tei) is extracted 
using the numerical model. The overall match of the 
magnitude of the cross section gives confidence in the 
numerical model. Error bars represent statistical varia- 
tion. In addition, there is systematic uncertainty due to 
uncertainty in the trapping potential of typically a factor 
of two, but this becomes more of an issue for lower sample 
temperatures and shallower traps, which are more chal- 
lenging to characterize and model due to the importance 
of gravity. The assumption of ergodicity may also be less 
valid at lower temperature because the elastic collision 
rate and rj become very small. 



VII. CONCLUSION 

In this paper, we have presented a model describing 
inelastic and elastic collision dynamics of trapped atoms 
that can treat traps lacking spatial symmetry and sam- 
ples with low r). The main assumptions are ergodicity 
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and a truncated Boltzmann velocity distribution. Statis- 
tical mechanical quantities calculated with this numeri- 
cal model were checked against analytic expressions for 
power-law traps. This also allowed us to discuss the range 
of validity of common assumptions made when describ- 
ing evaporation in power- law traps. The model was also 
used to describe ^*Sr in an asymmetric ODT with low 77, 
and coUisional parameters extracted from the data were 
found to agree well with theoretical predictions. 
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